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ABSTRACT 

Aims. The physical interpretation of spectro-interferometric data is strongly model-dependent. On one hand, models involving elab- 
orate radiative transfer solvers are too time consuming in general to perform an automatic fitting procedure and derive astrophysical 
quantities and their related errors. On the other hand, using simple geometrical models does not give sufficient insights into the physics 
of the object. We propose to stand in between these two extreme approaches by using a physical but still simple parameterised model 
for the object under consideration. Based on this philosophy, we developed a numerical tool optimised for mid-infrared (mid-IR) 
interferometry, the fast ray-tracing algorithm for circumstellar structures (FRACS) which can be used as a stand-alone model, or as 
an aid for a more advanced physical description or even for elaborating observation strategies. 

Methods. FRACS is based on the ray-tracing technique without scattering, but supplemented with the use of quadtree meshes and the 
full symmetries of the axisymmetrical problem to significantly decrease the necessary computing time to obtain e.g. monochromatic 
images and visibilities. We applied FRACS in a theoretical study of the dusty circumstellar environments (CSEs) of B[e] supergiants 
(sgB[e]) in order to determine which information (physical parameters) can be retrieved from present mid-IR interferometry (flux and 
visibility). 

Results. From a set of selected dusty CSE models typical of sgB[e] stars we show that together with the geometrical parameters 
(position angle, inclination, inner radius), the temperature structure (inner dust temperature and gradient) can be well constrained 
by the mid-IR data alone. Our results also indicate that the determination of the parameters characterising the CSE density structure 
is more challenging but, in some cases, upper limits as well as correlations on the parameters characterising the mass loss can be 
obtained. Good constraints for the sgB[e] central continuum emission (central star and inner gas emissions) can be obtained whenever 
its contribution to the total mid-IR flux is only as high as a few percents. Ray-tracing parameterised models such as FRACS are 
thus well adapted to prepare and/or interpret long wavelengths (from mid-IR to radio) observations at present (e.g. VLTI/MIDI) and 
near-future (e.g. VLTI/MATISSE, ALMA) interferometers. 

Key words. Methods: numerical, observational - Techniques: high angular resolution, interferometric - Stars: mass loss, emission- 
line, Be, massive, supergiants 



. | 1. Introduction medium density needs to be parameterised and it is not deter- 

, mined in a self-consistent way. For massive stars for instance, it 

r> ■ When dealing with optical/IR interferometric data, one needs to would be nece ssary to take into account non-LTE effects includ- 
ed ■ invoke a model for the understanding of the astrophysical ob- ing both gas and dust emission of the circumstellar material as 
ject under consideration. This is because of (1) the low coverage well as a m trea t me nt of radiation hydrodynamics. Fitting in- 
of the uv-plane and most of the time because of the lack of the terferometric data this way is as yet impossible because of corn- 
visibility phase, and (2) because our aim is to extract physical p U ting time limitations 
parameters from the data. This is particularly true for the Mid- 
Infrared Interferometric Instrument (MIDI, Leinert et al. 2003) Of course, solving at least the radiative transfer in a self- 
at the Very Large Telescope Interferometer (VLTI), on which consistent way is already very demanding for the computa- 
our considerations will be focused. Some pure geometrical in- tional resources. Consequently, model parameters cannot be de- 
formation can be recovered through a simple toy model such as termined in a fully automatic way and the model fitting process 
Gaussians (see e.g. Leinert et al. 2004; Domiciano de Souza must be carried out mostly by hand, or automatised by systemat- 
et al. 2007). ically exploring the parameter-space , the "chi-by-eye" approach 
However, this approach does not give any insights into the mentioned in Press et al. (1992). The followers of this approach 
physical nature of the object. One would dream of having a consider the "best fitting" model as their best attempt: a model 
fully consistent model to characterise the object under inspec- that is compatible with the data. It is admittedly not perfect, but 
tion. In many cases, if not all, a fully consistent model is out it is in most of the cases the best that can be done given the dif- 
of reach and one uses at least a consistent treatment of the ra- ficulty of the task. It is remarkable that a thorough x 1 analysis 
diative transfer. Models based for instance on the Monte Carlo of VLTI/MIDI data of the Herbig Ae star AB Aurigae has been 
method are very popular (see e.g. Ohnaka et al. 2006; Niccolini performed by di Folco et al. (2009) which remains to date one of 
& Alcolea 2006; Wolf et al. 1999) for this purpose. Still, the the most achieved studies of this kind. From the^- 2 analysis, for- 
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mal errors can be derived and at least the information concern- 
ing the constraints for the physical parameters can be quantified. 
Qualitative information about the correlation of parameters can 
be pointed out. 

The next step after the toy models for the physical char- 
acterisation of the astrophysical objects can be made from the 
pure geometrical model towards the self-consistency by includ- 
ing and parameterising the object emissivity in the analysis. 
For instance, Lachaume et al. (2007) and Malbet et al. (2005) 
use optically thick (i.e. emitting as black bodies) and infinitely 
thin discs to model the circumstellar environment of pre-main- 
sequence and B[e] stars. Of course this approach has some re- 
striction when modelling a disc: for instance it cannot handle 
nearly edge-on disc and an optically thin situation. 

We propose an intermediate approach: between the use of 
simple geometrical models and sophisticated radiative transfer 
solvers. Indeed, it is a step backwards from the "self-consistent" 
radiative transfer treatment, which is in most cases too advanced 
with regard to the information provided by the interferometric 
data. For this intermediary approach, we assume a prescribed 
and parameterised emissivity for the medium. Our purpose is to 
derive the physical parameters that characterise this emissivity. 
In the process, we compute intensity maps and most particularly 
visibility curves from the knowledge of the medium emissivity 
with a fast ray-tracing technique (a few seconds depending on 
map resolution), taking into account the particular symmetries 
of a disc configuration. Then, the model fitting process can be 
undertaken in an automatic way with standard methods (see e.g. 
Levenberg 1944; Marquardt 1963). The techniques we present 
are designed to be quite general and not tailored to any particular 
emissivity except for the assumed axisymmetry of the problem 
under consideration. 

Our purpose is twofold. On one hand - as already mentioned 
- we aim to estimate physical parameters and their errors char- 
acterising the circumstellar dusty medium under consideration 
with as few restrictive assumptions as possible; at least within 
the obvious limitations of the present model. On the other hand 
our purpose is to provide the user of a more detailed model, such 
as a Monte Carlo radiative transfer code, with a first characteri- 
sation of the circumstellar matter to start with. 

In Sect. 2 we describe the general framework of the proposed 
ray tracing technique. In particular how to derive the observable 
from the astrophysical object emissivity. In Sect. 3 we describe 
the numerical aspects that are specific to the present ray-tracing 
technique. In particular, the use of a quadtree mesh and the sym- 
metries that allow us to speed up the computation are detailed. 
In Sect. 4 we focus our attention on the circumstellar disc of 
B[e] stars and describe a parametric model of the circumstel- 
lar environment. In Sect. 5 we analyse artificial interferometric 
data generated both from the parametric model itself and from a 
Monte Carlo radiative transfer code (Niccolini & Alcolea 2006). 
Our purpose is not to fit any particular object, but to present our 
guideline to the following question: which physical information 
can we get from the data ? A discussion of our results and the 
conclusions of our work are given in Sects 6 and 7 respectively. 



2. The ray-tracing technique 

We describe here the FRACS algorithm, developed to study 
stars with CSEs from mid-IR interferometric observables 
(e.g. visibilities, fluxes, closure phases). Although FRACS could 
be extended to investigate any 3D CSE structures, we focus 
here on the particular case of axisymmetrical dusty CSEs. This 



4 




Fig. 1. Coordinate systems. The shaded ellipse represents a disc viewed 
by the observer. 



study is motivated by the typical data one can obtain from disc- 
like CSE observed with MIDI, the mid-IR 2-telescope beam- 
combiner instrument of ESO's VLTI (Leinert et al. 2003). 

2.1. intensity map 

Intensity maps of the object are the primary outputs of the model 
that we need to compute the visibilities and fluxes that are di- 
rectly compared to the observations. For this purpose, we in- 
tegrate the radiation transfer equation along a set of rays (ray- 
tracing technique) making use of the symmetries of the problem 
(see Sect. 3 for details). 

The unit vector along the line of sight is given by h = 
y sin i + z cos i, i being the inclination between the z-axis and 
the line of sight and x, y the unit vector along the x et j-axis of 
a cartesian system of coordinates (see Fig. 1), referred to as the 
"model system" below. The problem is assumed to be invariant 
by rotation around the z-axis. We define a fictitious image plane 
by giving two unit vectors Y = —y cos i + z sin i and X = -x. 
This particular choice is made making use of the axisymmetry 
of the problem. Note that for this particular coordinate system 
(X, Y) the disc position angle (whenever i + 0) is always defined 
as 90°. The actual image plane, with the Y' and X' axis corre- 
sponding respectively to North and East, is obtained by rotating 
the axis of our fictitious image plane by an angle PAj - |, where 
PAd is the position angle of the disc with respect to North. 

The dust thermal emissivity at wavelength A and position 
vector r is given by 

r U {r)=$Xr)B l {T{r)), (D 

where A^ bs (r) is the absorption coefficient and B,\(T(r) the Planck 
function at the medium temperature T(r) at r. is defined as 
n(r) CJj bs , where CJj bs is the absorption cross section and «(r) the 
number density of dust grains at r . 

We neglect the scattering of the radiation by dust grains, op- 
timising our approach to long wavelengths (from mid-IR to ra- 
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dio). This assumption simplifies the radiative transfer equation 
by removing the scattering term. 

We obtain the intensity map at position (X, Y) in the image 
plane (inclined by ;) and at wavelength A by integrating the trans- 
fer equation along the particular ray that passes through the con- 
sidered point of the image plane. Defining r s (X, Y, i) (simply r s 
for short) as the position vector along a ray, given in the model 
system of coordinates by 



r,(X,Y,i) 



-Y cos i + s sin i 
Y sin i + s cos ; 



(2) 



and by introducing the optical depth at wavelength A and posi- 
tion s along the ray by 



T X (X. 



\Y,i;s)= J 



(3) 



we obtain 



h(X,Y,i) 



■ / 



K f\r s )B A [T(r s )] e 



-T A (X,Y,i; s) 



ds , (4) 



where the extinction coefficient *^ xt (r) « K* s (r) because scatter- 
ing is neglected. 

We assume that the CSE is confined within a sphere of radius 

R out , s varies consequently from - ^Rl ut - R 2 to -^R^ - R 2 

(R 2 = X 2 + Y 2 ) in Eq. (4) and in the definition of a ray Eq. (2). 
This hypothesis can be relaxed without altering the present con- 
siderations and the domain of integration of Eq. (4) suitably cho- 
sen. 

If some radiation sources (e.g. black body spheres) are in- 
cluded in the analysis, an additional term must be added in 
Eq. (4) whenever a particular ray intersect a source. For a 
source with specific intensity Pj this additional term is given by 
P e -n(x,Y,i;s s ^ s (s) bgjjjg tne distance at which the ray given by X, 
Y and ;' (see Eq. 2) intersects the outermost (along the ray) source 
boundary. In that case the lower integration limit in Eq. (4), that 

is - yjRl ut - R 2 , must also be replaced by s (s) . 
2.2. Interferometric observables 

From the monochromatic intensity maps at wavelength A (Eq. 4) 
we obtain both the observed fluxes F A and visibilities V A for an 
object at distance d, 



FaH) 



CO oo 



hiX, Y, i) dXdY , 



(5) 



and 



3. Numerical considerations 

We seek to produce intensity maps within seconds 1 and we aim 
for our numerical method to be sufficiently general in order to 
deal with a large range of density and temperature structures. 
Given these two relatively tight constraints, the numerical inte- 
gration of Eq. (4) is not straightforward. 

For example we have tested that the 5 th order Runge-Kutta 
integrators of Press et al. (1992) with adaptive step-size (as dis- 
cussed in Steinacker et al. 2006) doest not suit our constraints. 
Indeed, the step adaption leads to difficulties if sharp edges 
(e.g. inner cavities) are present in the medium emissivity. 

3.1. Mesh generation 

Regarding the above mentioned constraints and the different nu- 
merical approaches tested, we found that Eq. (4) is more effi- 
ciently computed with an adaptive mesh based on a tree data 
structure (quadtrees/octrees). The mesh purpose is twofold: first, 
it must guide the computation of Eq. (4) and distribute the in- 
tegration points along the rays according to the variations of 
the medium emissivity; second, within the restriction of axis- 
symmetrical situations, the mesh must handle any kind of emis- 
sivity. Quad/octree meshes are extensively used in Monte Carlo 
radiative transfer codes (e.g. see Bianchi 2008; Niccolini & 
Alcolea 2006; Jonsson 2006; Wolf et al. 1999); the mesh gen- 
eration algorithm is thoroughly described in Kurosawa & Hillier 
(2001). 

The mesh we use is a cartesian quadtree. Cartesian refers 
here to the mesh type and not to the system of coordinates we 
use. Indeed, the mesh is implemented as a nested squared do- 
main (cells) in the p - \z\ plane (p = ^x 2 + y 2 ). The whole mesh 
is enclosed by the largest cell (the root cell in the tree hierarchy) 
of size R out in p and \z\. The underlying physical coordinate sys- 
tem is cylindrical (with z > 0) and the mesh cells correspond 
to a set of two (for z > and z < 0) tori, which are the actual 
physical volumes. 

The mesh generation algorithm consists in recursively divid- 
ing each cell in four child cells until the following conditions are 
simultaneously fulfilled for each cell in the mesh (see Kurosawa 
& Hillier 2001, for more details): 



Iff dh 
_Vf 

fff[<(r)] a J 3 ' 

v m 

fff [T(r)f dh 
Jff[T(r)fdh 



< n 



< n, 



and 



(7) 



(8) 



V A (B,PA)=—^ f fh(X,Y,i)e- 2} "-*^ cosW+ i sin(A) ] dXdY 

d 2 F A (i) J J 

— QO—CO 

(6) 

where V A is obtained for a given baseline specified by its pro- 
jected length B (on the sky, i.e. (X', Y') coordinates) and its polar 
angle PA from North to East (direction of the Y' axis). A and j 
represent, respectively, PAj - PA and V^T. 



where is the volume of cell V tot is the volume of the root cell 
and a, ft and rj are parameters controlling the mesh refinement. 

In the present work a and /3 have been fixed to 1, but higher 
values can be useful for some particular situations where the 
generated mesh must be tighter than the mesh generated directly 
from the /cj| bs and T variations. Typically, these situations show 
up for high optical depths (in this paper, optical depth values do 
not exceed ^ 1 at 10pm along the rays). The practical choice 
of a, p and 77 is obtained from a compromise between execu- 
tion speed and numerical accuracy of the Eq. (4) integration 



1 The actual computation time reached is less than 10 s for a 10 4 pixel 
map on an Intel T2400 1.83 GHz CPU. 
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Fig. 2. Quadtree mesh for a disc configuration. The disc parameters are 
those of model (b) described in Sect. 5.2 (see also Tables 2 and 4). The 
mesh refinement parameter 77 (Eqs. 7 and 8) has been set to the high 
value 1(T 3 in order to obtain a coarse mesh more easily represented. 



(e.g. typical values of 77 range from 10~ 5 to 10~ 4 ). When deal- 
ing with optically thick situations, a supplementary conditions 
can be added to Eqs. (7) and (8) in order to prescribe an upper 
limit to the cell optical depth. For instance, making use of the 
computation of the integral in Eq. (7), one can add the following 
criterion for cell £ (whose centre is (p^ , z%) and size A f ) 



1 



2np^A^ 



ff!< 



l (r)<Tr< A-riim , 



(9) 



where Arii m is the prescibed upper limit to the cell optical depth. 

For the moderate optical depths reached in this work, with 
values of 77 down to 1(T 5 and ATi; m set to 10~ 2 , the criteria of 
Eqs. (7) and (8) are the leading conditions to the mesh refine- 
ment. 

Figure 2 shows the mesh obtained in the particular case of a 
B[e] circumstellar disc (see Sect. 4) for models whose parame- 
ters are given in Table 2 (see caption for more details). 

The volume integrals in Eq. (7) and (8) are estimated by 
Monte Carlo integration. For a quantity f(r) and for the cell £ 
the integral Jff f(r) cPr is approximated by 
v t 



2n 



; / 

Zf— f P(—f 



2ttA| 



P f(p,z)dzdp * J] p k f(fi k ,z k ) , (10) 



k=\ 



where we made explicit use of the mesh coordinates and where 
(Pk, Zk) with k = 1 , • • ■ , N are chosen randomly and uniformly 
within the cell domain. 



3.2. Symmetries 

We can make use of the CSE symmetries to reduce the compu- 
tation domain of an intensity map from Eq. (4) to only a fourth 
of it and consequently reduce the computation time. 

Recalling the definition of a ray (Eq. 2), we have two notice- 
able identities for any disc physical quantity <D (e.g. /<^ bs , A^ xt , T, 
n, . ..) depending on r 



®(r s (X,Y,i)) 
*(r,(X,r,i)) 



®(r s (-X,Y,i)) and 
<D(r_,(X,-r,i)) , 



(11) 

(12) 



where Eq. (11) expresses the disc symmetry with respect to the 
y — z plane and Eq. (12) the point symmetry with respect to the 
origin of the model system of coordinates. 



3lane ana t,q. (IZ) the point symmetry with respect to the 
of the model system of coordinates, 
i iOm the above identities it is straightforward to deduce their 
counterpart for the intensity map 



h(X,Y,i) = h(-X,Y,i) and 

+£max 

- f Kf(r,,)ds' 

h(X,YJ) = h(X,-Y,i)e 



(13) 
(14) 



where s max = Y^ut ~~ Note tnat tne exponential factor in 
Eq. (14) has to be evaluated when computing h(X, Y, i) anyway; 
no extra effort is required to derive h(X, Y, i) from h(X, -Y, i) 
except for the multiplication of h(X, -Y, i) by this factor. 

3.3. Intensity map 

The fictitious image plane is split into a set of pixels whose po- 
sitions Xj and Y k are given by 



Xj = A x x; 



A F x \k + 



1 N 



1 _ 7Y\ 

2 ~ 2") ' 



(15) 
(16) 



where A^ = Ay is the pixel size in X and Y, and is the number 
of pixels in X and Y, and where 



< j,k< (N + 2) + 6, 



(17) 



where S = -1 for N even and 5 = otherwise and stands for 
the integer division. Taking into account the symmetries men- 
tioned in Sect. 3.2 only a fourth of the pixels need to be consid- 
ered. 

The evaluation of the integral in Eq. (4) is carried out for 
each pixel (Xj, Y k ) and along the ray r s (Xj, Y k , i). The intersec- 
tion points of the ray with the cell boundaries corresponds to a 
set of distances along the ray defined as 



so = 

Si = S1-1 + Asu for 1 < / < rt ce lls , 



(18) 

(19) 



where n ce ii s is the number of cells encountered along the ray, and 
Asi the distance crossed within the I th cell. 

We estimate numerically the optical depth t a (X, Y, i; s), de- 
fined in Eq. (3), via the midpoint rule quadrature by 



"cells -1 



n(X,Y,i;si)*if= J] Kf\r SM/2 )A Sl , 



(20) 



k=l 



where we defined sj+i/2 = Si + for / = 0, • • • , n ce ii s - 1. 
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The numerical estimate of I A (Xj, Y k , i) is obtained by 

"cells- 1 

h(X jt Y k ,i)* 2 < bS K,^^Ki, ; )) e " ?As '' (21) 

1=0 

The results for all j and k can then be obtained from the 
discrete counterpart of the symmetry relations (13) and (14) 



h(XN-i-j,Y k ,i) - h(Xj,Y k ,i) , 
h(Xj,Y N - k - U i) = h(Xj,Y k ,i)erW . 

3.4. Interferometric observables 



(22) 
(23) 



From the numerical estimate of I A (Xj, Y k , i) given above we ob- 
tain (similarly to Eqs. 5 and 6) the numerical fluxes and visibil- 
ities, which can be directly compared to the observed data. The 
numerical estimate of these quantities is again obtained through 
the mid-point rule. 

The numerical flux F A (i) is computed by 



1 



Fa(') * "T7 



V V 



d 2 

k=0 1=0 



h(X k ,Y h i) A X A Y . 



(24) 



The complex visibility is approximated numerically by 

. JV-1 N-l 

V * x Z Z M ' Y >> ^ AxAy ' (25) 

" rAW k=0 1=0 

where B = (B cos A, B sin A) and R kl = (X k , Y t ). 
3.5. Artificial data generation 

The procedure described below aims to mimic the observables 
of the VLTI/MIDI instrument: the flux F A (Eq. 24) and the mod- 
ulus of the visibility \V A \ (Eq. 25). The wavelengths and base- 
lines chosen for the artificial data generation correspond to ac- 
cessible values to VLTI/MIDI with the Unit Telescopes (UTs): 
Aj = 7,8,9, 10, 11,12, and 13//m(j = ,n A ;n A = 7), and 
(B k ,PA k ) as shown in Table 1 (k = 1, • • • ,n B ; n B = 18). These 
values amount to 126 points covering the uv-plane. 

For a given intensity map at Aj, F A . and \V Aj \ are taken as 
the expectation values of the simulated data. The observed flux 
Ff"* is then generated assuming a Gaussian noise with an RMS 
(root mean square) corresponding to 10 % relative error cr F ( j) = 
0.lxF Aj . 

The artificial observed visibility amplitudes |V° bs | are ob- 
tained as 



|V£"(B t ,PA*)| 



\V Aj (B k ,PA k )\ + AV k 



(26) 



where AV^ is a wavelength independent shift that mimics the 
error in the observed visibilities, introduced by the calibra- 
tion procedure commonly used in optical/IR interferometry. For 
each (B k ,PA k ), AV k is computed assuming a Gaussian noise 
with an RMS corresponding to 10% relative error (typical for 
VLTI/MIDI) o- v (k) = 0.1 x(\V A (B k , PA*)|>, where (\V A (B k , PA k )\) 
is the wavelength mean visibility modulus. 

3.6. Model fitting and error estimate 

We describe here the procedure adopted in order to simultane- 
ously fit observed fluxes and visibilities using FRACS models 



Table 1. Projected baselines. These values correspond to the baselines 
accessible from pairs of Unit Telescopes (UT) at ESO-VLTI. 



k 


B k [m] 


PA, [deg] 


1 


37.8 


61.7 


2 


41.3 


53.4 


3 


43.7 


44.8 


4 


46.2 


44.5 


5 


49.5 


37.5 


6 


51.9 


30.0 


7 


61.7 


134.6 


8 


62.0 


111.2 


9 


62.4 


122.5 


10 


81.3 


108.2 


11 


83.0 


52.2 


12 


86.3 


96.0 


13 


89.0 


84.4 


14 


89.9 


44.8 


15 


94.8 


36.7 


16 


113.6 


82.4 


17 


121.2 


73.6 


18 


126.4 


64.9 



defined by a given set of input parameters. This procedure is ap- 
plied to artificial data in the next sections. 

In order to quantify the discrepancy between the artificial ob- 
servations (|V^ S | and F°^ s ) and the visibilities and fluxes from a 

given model (|U^(B fc , PA^)| and F A .) we use the^ 2 like quantities 



"A n B 

in ~ ZZ 

7=1 k=\ 



|V*'(B it ,PA it )|-|V ily (B t ,PA*)r 



cr v (k) 



and 



j=\ k=\ 



F° DS - Fi 

A: 



o-fU) 



(27) 



(28) 



To take into account both the mid-IR flux and the visibilities 
on the same level in the fitting process, we minimise the follow- 
ing sum 



2 2,2 

X =X] V \ + Xf ■ 



(29) 



In the discussion below about the parameter and error deter- 
mination we use the reduced^ 2 defined by x\ = X 2 /(^- n B n A ~ 
"free) (for nfree f ree parameters). 

From a minimising algorithm the best-fit model parameters 
can then be found by determining the minimum y 2 : y 2 . . The 

j c /v r /v r.min 

"error" estimate is obtained from a thorough exploration of the 
parameter space volume, defined by a contour level ^ 2 min + A^- 2 , 

where A^ 2 has been chosen equal to 1 . This volume can be inter- 
preted as a confidence region. The quantity defined in Eq. (29) 
is a weighted sum of x 2 variables whose cumulative distribu- 
tion function can be approximated by a gamma distribution (see 
Feiveson & Delaney 1968) with the same mean and variance. 
It is then possible to obtain a rough estimate of the confidence 
level associated with the Ax 2 = 1 confidence region given ap- 
proximately by — 2 o~. 

The size of the confidence region is determined by consid- 
ering all possible pairs of parameters for a given fitted model 
and computing x 2 maps for each. The procedure to estimate the 
errors can be summarised as follows: 

- For a given x 2 map, i.e. for a given couple of parameters 
among the «f ree x («f ree - l)/2 possibilities, we identify the 
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Mid-IR flux (model a) 




Mid-IR flux (model b) 




Fig. 3. Simulated VLTI/MIDI mid-IR object flux. The numerically generated data (left model a, right model b) are shown as circles with error 
bars. The values of the fluxes for the best-fit models are represented as crosses. 



region bounded by = 1 around the minimum of this 
particular map. 

- The boundaries of the projection of these regions on each 
of the two parameter axis considered are recorded for each 
map. 

- The final errors on a given parameter are taken as the highest 
boundary values of the projected regions over all maps. 



4. Astronomical test case: sgB[e] stars 

In the following sections we apply FRACS to a theoretical in- 
terferometric study of dusty CSE of B[e] supergiants (sgB[e] in 
the nomenclature of Lamers et al. 1998). However, we empha- 
sise that FRACS is in no way restricted to this particular class of 
objects. 

sgB[e] stars reveal in particular a strong near- or mid-IR 
excess caused by hot dust emission. There is evidence (e.g. 
Zickgraf et al. 1985) that the stellar environment, and in par- 
ticular dust, could be confined within a circumstellar disc. Our 
purpose is to characterise this class of objects and derive not 
only geometrical parameters (e.g. inner dust radius, disc posi- 
tion angle and inclination) but also physical parameters such as 
temperature gradients, dust formation region, material density, 

The physical description of the CSE chosen for our study is 
the wind model with equatorial density enhancement. This is a 
classical CSE model commonly adopted for sgB[e] (e.g. Porter 
2003). 

In order to compute the model intensity maps we need to pa- 
rameterise the emissivity of the disc. Consistently with FRACS 
assumptions, we consider only dust thermal emission without 
scattering by dust grains and the gas contribution to the medium 
emissivity. In the rest of this section we characterise the emis- 
sivity by describing the dust density law, the absorption cross 
section, and the temperature structure of the CSE. 



4. 1 . Mass loss and dust density 

Dust is confined between the inner and outer radius R m and R out 
respectively. We assume a stationary and radial mass loss; phys- 
ical quantities will consequently depend only on the radial com- 
ponent r and the co-latitude 9. The disc symmetry axis coincides 
with the z axis of the model cartesian system of coordinates. The 
mass loss rate and velocity parametrisations are simplifications 
of the one adopted by Carciofi et al. (2010), and we refer the 
reader to their work for a complete description (see also Stee 
et al. 1995, for a similar description). 

The mass loss rate per unit solid angle, at co-latitude 9, is 
parameterised as follows 



^(0) = ^(0X1+^1 sin m (0)) . 
all ail 



(30) 



with the help of two dimensionless parameters A\ and m. 

Even though our computations make no explicit use of the 
radial velocity field v r (ff) (assumed to have reached the termi- 
nal velocity v 00 (6') in the region under considerations, i.e. v r {8) « 
Vco(0)), the dust density depends on v r {9) parameterised in a sim- 
ilar fashion 

v r (ff) = v,.(0) (1 +A 2 sin'" 6) , (31) 

where we have introduced the supplementary dimensionless pa- 
rameters A 2 . From Eqs. (30) and (31) we see that A\ and A2 are 
the relative differences of the values of ^(9) and v r (9) at the 
equator and the pole (relatively to the pole). 

From the mass continuity equation one obtains the number 
density of dust grains 



«(r,0) = n in (— j — 



1 +A 2 l+Ai (sin0)" 



Ai i+a 2 (sin ey 



(32) 



where n m is the dust grain number density at Ri n in the disc equa- 
torial plane. In Eq. 32, the parameter m controls how fast the 
density drops from the equator to the pole, defining an equato- 
rial density enhancement (disc-like structure). 
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Consistent with the accepted conditions for dust formation 
(Carciofi et al. 2010; Porter 2003) we assume that the dust can 
survive only in the denser parts of the disc. We thus define a 
dusty disc opening angle AO^ determined by the latitudes for 
which the mass loss rate has dropped to half of its equatorial 
value: 



A#d = 2 arccos 



2 A x 



(33) 



To summarise, the dust grains only exist (i.e., n(r, 9) + 0) in 
the regions bounded by R m <r< R om and by < g < 

4.2. Dust opacities 

The absorption cross section C'f s for the dust grains is ob- 
tained from the Mie (1908) theory. The Mie absorption cross 
sections are computed from the optical indices of astronomical 
silicate (Draine & Lee 1984). Note that since scattering is ne- 
glected, CJj bs « CJj xt , with C™' being the extinction cross section. 

For a power-law size distribution function according 
to Mathis et al. (1977) the mean cross sections (e.g. for C'f s ) 
are given by 



jf a^q bs (fl)d/fl 



abs min 



a max 



(34) 



/ a^da 



Table 2. Model parameters. This table lists the parameters of 12 differ- 
ent models. The parameter values that change (A6> d , n m and i) from one 
model to the other have been inclosed in a box and separated by a slash. 
The values of A6 A = 10°/60° given below correspond via Eq. (33) to 
m = 183.56/4.86 respectively. 



Parameters 



Values 



Unit 



A, 
A 2 
M A 

Rm 

Rout 

n m 
T m 
7 

P M 
a 

PA d 
i 

Hrnax 



150 
-0.8 



10/60 



60 
30 
3000 



0.015/0.15 



1500 
0.75 
6500 

3 

125 



20/50/90 



0.5 
50 
-3.5 



deg 
R Q 
R s 
R s 
m- 3 
K 

Wm^/iirr 1 str~ 

deg 
deg 
Hm 



to lie between a — 4 (pure black body) and a ^ 2.6 (free-free 
emission) for an electron density proportional to r~ 2 (Panagia & 
Felli 1975; Felli & Panagia 1981). 



where a m [ n and a max are the minimum and maximum radii for 
the dust grains under consideration and /3 is the exponent of the 
power-law. The computation of the cross section in Eq. (34) was 
performed with the help of the Wiscombe (1980) algorithm. 

4.3. Temperature structure 

The dust temperature is assumed to be unique (i.e. independent 
of grain size) and described by a power-law 



T(r) = T h 



(35) 



where 7^ is the temperature at the disc inner radius R{ n . We note 
that y is not necessarily a free parameter because in the optically 
thin regime (large wavelength and radius) the temperature goes 

2 

as T(r) oc r~~s with 6-1 (see Lamers & Cassinelli 1999). 
4.4. Central continuum emission 

The continuum emission from the central regions is composed 
by the emission from the star and from the close ionised gas 
(free-free and free-bound emission). This central source emis- 
sion is confined to a small region of radius /? s («: R[ n ), which is 
unresolved (angular sizes of a few milliarcseconds) by mid-IR 
interferometers. Thus, in our modelling R s is simply a scaling 
factor of the problem fixed to a typical radius value for massive 
stars. The specific intensity (in Wm^yUtrT 1 str -1 ) of this central 
source is parameterised as follows 



1 



(36) 



where I* is the specific intensity at a reference wavelength A Q 
(=10 //m in the following), and a gives the spectral dependence 
of the continuum radiation. In the mid-IR its value is expected 



5. Study of the tested models 

Following the description in the last sections we describe 
here the chosen sgB[e] model parameters used to simulate 
VLTI/MIDI observations (visibilities and fluxes) and the cor- 
responding analysis, i.e. model fitting, using FRACS. The list 
of chosen parameters is summarised in Table 2. Two types of 
numerical tests are presented. Firstly, synthetic mid-IR interfer- 
ometric data are generated from FRACS itself. In that way, it 
is possible to determine what information the mid-IR interfero- 
metric data contain under the optimistic assumption that we do 
have the true model. Secondly, this study is supplemented by the 
comparison of FRACS to a Monte Carlo radiative transfer com- 
putation. This confirms that FRACS can indeed mimic, under 
appropriate conditions, the results of a more sophisticated code 
as seen from the mid-IR interferometric eye. 



5.1. Parameter description 



1 kpc, 



The distance to the simulated object has been fixed to d ■ 
which is a typical distance for Galactic sgB[e]. 

The inner radius R m = 30 R s = 1800 Rq value was cho- 
sen by considering the location of the hottest dust grains (see 
Lamers & Cassinelli 1999) with a condensation temperature of 
1500 K assumed to be the T m value. The value of R oyll can- 
not be determined from the mid-IR data and has been fixed to 
3000 R s = 1.8 x 10 5 Rq. The temperature gradient y was fixed 
to 0.75 according to Porter (2003). PA d was fixed to 125 °. 

The central source emission is supposed to have a radius 
R s = 60 Rq. We recall that the central region is unresolved by 
the interferometer and that its radiation describes both the stellar 
and inner gas contribution to the continuum mid-IR emission. 
The specific intensity of this central source I s , has been chosen 

to be 6500 Wm 4 /im 4 str -1 . If the central source was a pure 
blackbody this value would correspond to the 10;um emission 
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of a blackbody with an effective temperature around 8000 K. 
However, this central emission is not a pure blackbody, and we 
adopt the spectral dependence of the central source emission to 
be a = 3, which is a compromise between a = 4 for a pure 
blackbody and a value of 2.6 for free-free emission (Panagia 
& Felli 1975; Felli & Panagia 1981). 

Spectroscopic observations of Ha and forbidden line emis- 
sions from B[e] CSE (Zickgraf 2003) reveal that typical values 
for A 2 are expected to range from -0.95 to -0.75. We adopt 
the value -0.8 in our models. According to Lamers & Waters 
(1987), the values of A\ range from 10 2 to 10 4 in most cases 
(though values as low as 10 are not excluded). With this high 
value of A] the factor 1 + Ai(sin6>) 1/m in Eq. (32) of n(r, 6) is 
approximatively given by A 1 (sin Q) l l m for all pertinent values of 
6, i.e. those close to n/2 within the disc. This leads to an evi- 
dent degeneracy in n in x A\ in n(r, 6): we are only sensible to the 
product of the two parameters as a scaling factor for the density. 
Therefore, the value of A\ is assumed to be fixed to 150. 

To define the dust opacities the chosen value for /3 is that 
of Mathis et al. (1977), i.e./? = -3.5. Because some sgB[e] 
show weak 9.7 //m silicate features in their spectrum (e.g. Porter 
2003; Domiciano de Souza et al. 2007) we chose to use large 
grains in our test models: a m ; n = 0.5 jum and a max = 50//m. 
However, with this particular choice of large grains, the average 
albedo from 7 to 13 [im is 6.4 %, with the highest value reached 
at 7 //m. We have checked with a Monte Carlo (MC) simulations 
(see Sect. 5.3) that the effect of scattering on our primary ob- 
servables, visibilities, and fluxes is indeed negligeable by com- 
paring the results obtained by switching the scattering process 
off and on 1 . The mean relative differences are 3.5 % and 3.0 % 
for the visibilities and the fluxes respectively. These values must 
be compared to the effect of random noise in the MC simulation, 
estimated to be of the same order and to experimental errors, 
typically ~ 10 % for the visibilities and fluxes. We underline that 
whenever the albedo can be neglected, it is theoretically safe to 
compute visibilities and fluxes from the consideration presented 
in Sect. 2, in any other situations the effect of scattering on the 
observable must be carefully tested. 

The parameters n m , m and i were set to different values 
defining 12 test models to be analysed from their correspond- 
ing simulated data. Two n m values (0.015 irT 3 and 0.15 trT 3 ) 
have been chosen in order to have an approximate disc-dust 
optical depth in the equatorial plane (from R m to R out ) close 
to =! 0.1 and ^ 1 in the wavelength range considered (from 
7/im to 13//m). These values corresponds to a mass loss rate 
of M =s 2.5 x 10~ 7 ' "~ 6 Mq yr _1 . Two m values were chosen cor- 
responding to a wide and a narrow opening angle, i.e. Aflj = 10 
and A#d = 60°. Three inclinations i were tested (20°, 50°, and 
90°) corresponding to discs seen close to pole-on, intermediate 
inclination, and equator-on. These values of n m , m, i, together 
with the parameters fixed above, define 12 test models that will 
be studied below. 

From these 12 test models we have generated 12 sets of arti- 
ficial VLTI/MIDI observations (visibilities and fluxes) following 
the procedure described in Sect. 3.5. We do not aim to present an 
exhaustive revue of all types of sgB[e] CSE. Rather, we focussed 
on the analysis of the parameter constraints one can hope to ob- 
tain from present and near-future mid-IR spectro-interferometry. 
The quantitative estimate of these constraints is derived from a 
systematic analysis of the^f 2 variations with the parameters. 



2 The computation have been done for model b described in Sect. 5.2, 
the baselines listed in Table 1 and the wavelengths under consideration 
from 7 to 13jum. 



In our model fitting and xl analysis we concentrate on 
10 free parameters («f ree = 10) that can be set into four different 
groups: 

- The geometrical parameters : PA d , i and R in , 

- the parameters related to the central source : /? and or, 

- those describing the temperature structure : T m and y, 

- and the number density of dust grains : A2, A#a (or equiva- 
lent^ m) and n; n . 

The remaining parameters of the model are in general 
loosely constrained by mid-IR interferometric observations so 
that we kept them fixed to the values described above. 

5.2. Model fitting andx\ analysis of the 12 test models 

We describe here the data analysis procedure adopted to study 
our 12 test models. The results of our analysis are summarised 
in Tables 3 and 4, and their physical interpretation is presented 
in Sect. 6. 

As a first step we chose 2 of the 12 models, hereafter called 
models (a) and (b), to be exhaustively studied from a complete 
model fitting procedure. As an example we show the simulated 
observed mid-IR fluxes and visibilities for model (a) and (b) in 
Figs. 3, 4, and 5. The parameters of models (a) and (b) are those 
of Table 2 with A6» d = 60 °, i = 50 and n in fixed to the val- 
ues 0.015 trr 3 and 0.15 irT 3 respectively. These two models are 
those presenting some of the best constrained model parameters 
for the dust CSE. On the other hand, the contribution of the cen- 
tral regions to the total flux and visibilities is quite different in 
models (a) and (b) (see discussion in Sect. 6). 

The study of models (a) and (b) have thus been performed 
as for real interferometric observations. The best-fit values of 
the parameters have been obtained by the Levenberg-Marquardt 
algorithm with a stopping criterion corresponding to a relative 
decrease in^ 2 of 10~ 3 . 

The errors on each model parameter have been obtained fol- 
lowing the methods described in Sect. 3.6. The \\ maps have 
been computed with a resolution of 21 x 21 around the best- 
fit values of the parameters. The map sizes have been adjusted 
in order to enclose the A^ 2 = 1 contour. This adjustment was 
performed until an upper limit for the map size of 100 % of the 
best-fit parameter values was reached. This amounts to the com- 
putation of 3.969 x 10 4 different models. The results, namely the 
mean relative error up to 100 %, for these two particular models 
are summarised in Table 3. 

The other ten models (numbered from 1 to 10 in Table 4) 
have been used in order to get some quantitative (but limited) 
information about how the uncertainties of the fitted parame- 
ters evolve as a function of three disc characteristics: its optical 
depth (t by means of n,„ parameter), its inclination (;) and its 
opening angle (A6<s, controlled by rti). To perform this study we 
have decided to limit the exploration of the space parameter in 
a relative range of 25 % on both sides from the model param- 
eters. In order to reduce the computation time, the maps were 
not generated around the best-fit parameters which would have 
required to compute several thousands models more but around 
the true parameters themselves. This procedure has the supple- 
mentary advantage that we do not rely on any specific minimi- 
sation algorithm. We checked that estimating the best-fit param- 
eters from the true ones is reasonable within a few percents us- 
ing the Levenberg-Marquardt algorithm with a stopping criterion 
corresponding to a relative decrease in^- 2 of 10~ 3 . The resolution 
of the^ 2 maps have been reduced to 15 x 15. The total number 
of models to be computed is as large as 1.0125 X 10 5 . 
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Table 3. Relative errors (given in %) on the parameters for models (a) et (b) (see text for description of models). For each of the 10 free parameters 
considered in the analysis, the values of the relative error corresponding to the 12 different models are given. Indeed, these relative errors are 
"mean values" for the errors because the error bars are not symmetric with respect to the best-fit values. The parenthesis around the relative error 
of A 2 recall that this parameter is bounded. 



Models \Parameters 


A 2 


m 




Rin 


"m 


T m 


7 




a 


PA d 


i 


(a) 


(53) 


> 


100 


1.9 


46 


15 


7.1 


27 


13 


4.3 


6.2 


(b) 


(59) 


> 


100 


4.5 


100 


20 


12 


> 100 


96 


9.4 


7.2 



Table 4. Constraints on the model parameters. For the 12 models considered here (differing in their value of t, A9 a and i), numbered from 1 to 10 
(except for model a and b), we classified the parameters into 3 different relative error ranges : below 10 %, between 10 and 25 % and above 25 %. 
Because A 2 , m, n- m , I\ q are determined for all the models with an error greater than 25 % they have been discarded from the table for the sake of 
clarity . 





Models \ parameters 






Constraints [%] 




T 


A(? d [deg] i 


[deg] 


< 10 


10 -> 25 


> 25 


1 


0.1 


60 


20 


Rm,7 


7m, « 


PA d , i 


2 


1. 


60 


20 


Kin 


Tin, 7 


a, PA d , i 


3 


0.1 


10 


20 




Rin 


T^, 7, <*, PAd, i 


4 


1. 


10 


20 


Rin, 7 


T m 


a, PA d , i 


(a) 


0.1 


60 


50 


Rin, 7, PAd, i 






(b) 


1. 


60 


50 


Rin, PA d , i 


7m, r 


a 


5 


0.1 


10 


50 






T m , 7, a, PA d , 


6 


1. 


10 


50 


Rin, y, PA d , ( 




a 


7 


0.1 


60 


90 


Rin, 7, PAd 


Tin, i, a- 




8 


1. 


60 


90 


Rin, PA d 


7m, 7, i 


a 


9 


0.1 


10 


90 


PAi 


Rin 


T^, 7, «, ; ' 


10 


1. 


10 


90 


Rin, 7, PA d , i 


T m 


a 



Table 5. best-fitting FRACS parameters from artificial data generated with a Monte Carlo code. The column "true values" refers to the MC input 
parameters, except for T- m , y, and y' which are determined from the results of the MC simulation. The columns "two power-law" and "one power- 
law" list the best-fit parameters obtained with FRACS assuming two and one power-law for the temperature respectively. The xlmin vames are 
respectively 0.73 and 0.79 for two and one power-law. 



Parameters 


Units 


true values 


two power-law 


one power-law 


A 2 




-0.8 


-0.791 


-0.782 


m 




4.86 


5.59 


4.74 


Rin 


R, 


30 


29.8 


29.9 


"in 


nr 3 


0.15 


0.189 


0.169 


T 

1 in 


K 


1150° 


1090 


1070 


ylY 




0.725/0.478° 


0.719/0.613 


0.676 




Rin 


5.24° 


2.87 




PA d 


Wm^yunT 1 str~' 


6.62 x 10 3 


6.48 x 10 3 


5.04 x 10 3 


deg 


125 


125 


124 


i 


deg 


50.0 


50.6 


50.2 



" These values are not prescribed parameters, but are determined from the results of the Monte Carlo simulation. The values reported here are 
best-fit parameters of the mean disc temperature (see text for more details). 



5.3. Comparison with a Monte Carlo simulation 



We generated synthetic data with the help of a Monte Carlo 
(MC) radiative transfer code (Niccolini & Alcolea 2006) for 
model b (see above) for the seven wavelengths considered in the 
problem and the baselines of Table 1 . Again, the adopted pro- 
cedure to generate the mid-IR interferometric data follows the 
considerations of Sect. 3.5. In the MC code, the source of pho- 
tons is described by a blackbody sphere of radius R s = 60 Rq 
and an effective temperature of T e s = 8000 K. The temperature 
of the CSE is not prescribed but computed from the Lucy (1999) 
mean intensity estimator. This choice of T e g gives at the inner 
radius of model b a dust temperature of ^ 1 150 K lower than the 
sublimation temperature. In this way, we can test if in the fitting 
process using FRACS, a spurious effect might not lead the min- 



imisation algorithm to reach the upper limit for T m of 1500 K, 
corresponding to the adopted dust sublimation temperature. 



We obtained the best fitting parameters for the CSE model 
described in Sect. 4 with FRACS. For a comparison with the MC 
code, a has been set and fixed to 4 corresponding to the value of 
a blackbody. Depending of the disc optical depth, the tempera- 
ture structure may show two separate regimes corresponding to 
(1) the inner regions with the strongest temperature gradient, op- 
tically thick to the stellar radiation and (2) the outer regions opti- 
cally thin to the disc radiation with a flatter temperature gradient. 
In order to determine if mid-IR interferometric data are sensitive 
to two temperature regimes, we tested the effect of two parame- 
terisations of the temperature structure: the unique power-law of 
Eq. (35) and a generalisation to two power-laws with a transition 
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radius, Rj, and a second exponent y' 
for r > Rj. 

The best-fitting parameters for both parameterisations are 
shown in Table 5. The images of the disc at 10 /im generated 
with the MC code and their corresponding FRACS counterpart 
(best-fitting model) are shown in Fig. 6 for comparison. 



6. Discussion 

We first discuss the uncertainties in the parameters derived for 
the 12 models studied in Sect. 5.2. For each model we divided 
the parameters into three groups associated to a given level of 
constraints expressed by the relative errors: below 10 %, between 
10 % and 25 %, and above 25 %. This information is summarised 
in Table 4. The exact relative errors for the two models studied 
in detail (models a and b) are shown in Table 3. Then, we anal- 
yse the results of Sect. 5.3 obtained from the best-fit of the data 
simulated with the MC radiative transfer code. 

6.1. Central source 

Table 4 shows that the central source parameters (I s , and a) 
can only be constrained with relative uncertainties > 10% for 
all test models. A deeper and more quantitative investigation of 
these parameters can be obtained from models (a) and (b). From 
Table 3 it can be seen that and a are much better constrained 
for model (a) (27 % and 13 %, resp.) than for model (b)(relative 
errors 100%). 

The key quantity for a good constraint for the central source 
parameters (/? and a) is simply the relative contribution of the 
flux of the central source to the total flux of the object (source 
and disc). Indeed, the models in Table 3 only differ by this rela- 
tive flux contribution of 5.3 % in model (a) (r = 0.1), while it is 
only 0.7 % in model (b) (t = 1). 

Our analysis thus shows that interferometric data can con- 
strain and a with a relative precision of ^ 15 % - 30 % even 
when the central source contributes to (only) a few percent of 
the total mid-IR flux. 

6.2. Geometrical parameters 

The parameters PAd, i, and Ri n are those usually estimated from 
simple geometrical models (e.g. ellipses, Gaussians). However, 
their determination from geometrical models is quite limited, in 
particular for i, for which only an estimate can be derived from 
the axis-ratio of an ellipse, for example. In addition, the esti- 
mate of ; from a simple analytical model such as a flat ellipse is 
only valid for configurations far from the equator (intermediate 
to low /). The use of a more physical and geometrically consis- 
tent model such as FRACS allows us to relax this constraint and 
makes the determination of i possible for all viewing configura- 
tions. 

As expected, PAd and ; are better determined if the inclina- 
tion of the disc with respect to the line of sight is away from 
pole-on (high ;'). In Fig. 8 we can clearly see this behaviour from 
the^ maps involving PAd and ;. Moreover, the uncertainties on 
PAd and ; do not seem to be strongly dependent on t (equiva- 
lently n m ) and A#d (equivalently m) for all models. 



The inner dust radius R m is not strongly dependent on any 
parameter (t, A#d or i), being very well constrained (better than 
10 %) for most tested models. 

6.3. Temperature 

The parameters related to the temperature structure of the CSE, 
r in and y, are well constrained in most models, with relative er- 
rors below 20 % and 12 % for both models (a) and (b). Indeed, 
y has a strong impact on the IR emission across the disc, and 
consequently this parameter has a direct influence on the visibil- 
ities (see Figs. 4, 5). T; n has a lower influence, compared to y, 
on the shape of the monochromatic image (radial dependence of 
intensity) and can be mainly considered as a scaling factor to it. 
On the other hand, the mid-IR flux imposes stronger constraints 
on Ti n . From Table 4 we see that the CSE's temperature structure 
is not highly dependent on t (n ln ) and A#d for tested models. 

6.4. Number density of dust grains 

The parameters related to the density law, that is to say m, n m 
and A2, seems to be rather poorly constrained from the mid-IR 
data alone. From the results of model (a) and (b) corresponding 
to an intermediary inclination i = 50 °, we found that only n m is 
constrained somewhat moderately with a mean relative error of 
46 %. For m and A2, according to the results of Table 3 it seems 
that nevertheless, upper limits to their values can be determined. 
Note that because A 2 is bounded (-1 < A2 < 0) the mean rel- 
ative errors, 53 % and 59 % for model (a) and (b) respectively, 
correspond approximatively to the limit values of A2, which is 
consequently not constrained. 

Table 4 confirms this trend for m, n m and A2 at least for the 
situations explored via the models presented here. From all maps 
computed within +25 % of the true value, we always found that 
the mean relative error to these parameters is larger than 25 % 
with no hint that it could be close to these limits. 

From Fig. 9, comparing the;^ maps for all pairs of n m , m and 
A2, we can see that the tyft contours get sharper around the min- 
imum value for model (a) (corresponding to lower optical depths 
along the line of sight) than for model (b). Indeed, the constraints 
on «i n and m are improved for lower optical depths, or equiv- 
alently for lower disc masses. Indeed, when the disc mass (or 
optical depth) decreases, the flux (mid-IR flux, intensity maps) 
emitted by the disc reflects the mass of the disc, while for high 
optical depths we only probe the regions of the disc very close 
to the projected surface revealed to the observer. A2, however, 
is unaffected by the change in disc mass and remains undeter- 
mined anyway. From Fig. 9 it can be seen that n in , m and A 2 are 
strongly correlated. This is expected from the expression of the 
density (see Eq. 32) depending on these parameters. However, 
this dependence and the final correlation between these param- 
eters are related through the computation of the visibilities and 
the mid-IR flux, as well as the comparison to the data and is, 
therefore, not straightforward. 

To improve the situation concerning n m , m and A2, the mid- 
IR data can be supplemented by other types of observations such 
as for instance spectroscopic data, from which one can better de- 
termine A2 (e.g. see Chesneau et al. 2005). We tested the effect 
of fixing the value of A2, or equivalently of assuming that A2 is 
fully determined, in the process of estimating the errors of the 
other parameters. For model (a), the relative errors on n m and 
m go down to 33 % and 71 % respectively while for model (b), 
«in and m are determined with an accuracy reaching 95 % and 
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78 % respectively. The precision to which other parameters are 
determined is not affected by the determination of A%. 

We also tested the influence of the determination of «i n , m 
and A2 on other parameters by fixing their values and estimatng 
the relative errors on the remaining parameters for model (a) 
and (b). Only I\ and T± a are more strongly affected by the deter- 
mination of «i n , m and A2: for model (a) (resp. model b) I s lg gets 
determined down to 19 % (resp. 80 %) and T m down to 9 % (resp. 
18 %). The influence is stronger with lower disc mass (model a 
compared to model b). This effect can be explained because if 
we have a good determination of the disc mass because we know 
n m , m, and A2, the determination of the parameters that scale the 
source and disc fluxes is improved accordingly for the visibility 
and the mid-IR flux. 

n m , m and A2 shape the density structure of the circumstellar 
medium. Though they are not well constrained, they certainly 
have a strong influence on the temperature structure, which in 
turn is very well constrained. For the particular case of sgB[e] 
circumstellar discs, a natural evolution of FRACS is to include 
the direct heating of the medium by the central source of radia- 
tion assuming that the disc is optically thin to its own radiation. 
The temperature structure would not be parameterised, and its 
good determination would certainly put better constraints on n; n , 
m and A2, while keeping an affordable computation time for the 
model-fitting procedure. This will be the purpose of a subsequent 
work. 

Finally, one can derive a ranking of the parameter constraints 
according to two criteria: first the parameter must be constrained 
within the prescribed limits (100 % for model a and b and 25 % 
for model 1 to 10) and second the mean relative error must be as 
low as possible. The best-fitted parameters, most of the time ac- 
cording to these criteraria are by decreasing order of best deter- 
mination: R[ n , PAd, y, Ti n , i, a, P., m and A2. This tendency 
can be seen in Table 4. 

6.5. best-fit to the MC simulation 

The Xrmm values obtained for the two types of temperature 
parametrisations (one and two power-laws) are quite similar: 
0.79 and 0.73 respectively. Regarding the data, both tempera- 
ture parametrisations are indeed acceptable. In addition, these 
results show that we can actually obtain very good fits from 
data sets based on more physically consistent scenarios. A com- 
plete error analysis and study of the parameter determination has 
been presented in the previous sections for data generated from 
FRACS and will not be repeated here. In particular, parameter 
confidence intervals, from which errors were derived, have al- 
ready been estimated. Here, we will instead focus on the true er- 
rors, i.e. the differences between the true model parameters and 
the best-fitting values for the parameters (see Table 5). The two 
types of errors must not be confounded. The true errors reflect 
the capability of FRACS to mimic the mid-IR interferometric 
data regarding the information it provides. Of course, with the 
sparse uv-plane coverage inherent to this kind of data as well as 
the experimental noise, one should not expect a full agreement 
of the fitted and the true parameters: they are indeed different. 

From Table 5, we see that the geometrical parameters, PAd, i 
and Rin, can be almost exactly recovered as expected. The source 
specific intensity /* and the parameters related to the density, 
A2, «in and m can be recovered fairly well and have best-fitting 
values close to the true parameters. 

The values of T\ n , R-y, y and y' reported as "true" in Table 5 
are indeed the values of a fit to the average (over the co-latitude 




Fig. 7. Temperature of the CSE. The solid line represents the best fit 
(last column in Table. 5) with a unique power-law, the dashed line the 
best-fit with two power-laws (fourth column in Table 5) and the dot- 
dashed line the MC results. The shaded region represents the possible 
domain for a unique power-law by taking into account the errors esti- 
mated in Table 3. 



for a given r) computed temperature in the disc. The true rela- 
tive differences for T m do not exceed 7 % independently of the 
adopted parameterisation of the temperature (one or two power- 
laws). The best-fitting values of y, the inner temperature gradi- 
ent, obtained with FRACS are very close to the true values with 
two and one power-law with true relative error of 1 % and 7 % 
respectively. This already suggests that the mid-IR data provide 
information on the inner and hottest region of the CSE, in par- 
ticular on the inner temperature gradient y. 

Fitting the temperature computed with the MC code with 
a simple power-law, we obtain y 0.64. This value is close 
to those of the best fitting models, especially with a unique 
power-law (6 % relative difference). For comparison, the actual 
mean temperature gradient as derived from the MC simulation 
is =! 0.60. For this particular data set, the values of y' and Rj 
recovered by FRACS differ by 28 % and 45 % respectively from 
the actual values. This again confirms the sensibility of the inter- 
ferometric data to the temperature structure mostly in the inner 
(r S= Rj) regions of the disc. The best-fitting models (fourth and 
last columns in Table 5) as well as the MC results are shown in 
Fig. 7. Regarding the errors (estimated from the results given in 
Table 3) shown as a shaded area, we can see that both tempera- 
ture parameterisations are essentially the same and show a better 
agreement with the MC results in the inner than in the outer re- 
gions of the disc. 

We considered a "truncated" model with two power-laws 
(with parameter values listed in the third column of Table 5) in 
which the CSEs emission for r > Rj, the "outer" regions, has 
been set to 0. We then compared the visibilities and the fluxes 
of this truncated model to the same model including the outer 
region emission. We obtained relative differences, averaged over 
all considered wavelengths and baselines (Table 1), of 18 % and 
17 % respectively. These relative differences are larger, but are 
still close to the noise level. For this reason, one cannot expect 
to obtain much information on the outer temperature gradient y', 
at least for the particular configuration we considered. 
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7. Conclusion 

We proposed and described here a new numerical tool to inter- 
pret mid-IR interferometric data. Even though we focussed on 
the special case of circumstellar disc observations, the numeri- 
cal techniques have been developed with the aim to be as general 
as possible. The methods we employ rely on both parameterised 
physical models and the ray-tracing technique. The need for such 
a tool is evident because the nature of interferometric data im- 
poses an interpretation through a model of the object to obtain 
any kind of information. On one hand, Monte-Carlo radiative 
transfer methods require too much computation time to asso- 
ciate the model-fitting to an automatic minimum search method. 
On the other hand, purely geometrical function fitting (such as 
ellipses or Gaussians) are too simple to envisage to obtain phys- 
ical constraints on the observed disc. Hence, a tool like FRACS 
fills a blank in the model fitting approach for mid-IR interfero- 
metric data interpretation. The main advantages of FRACS are 
its speed and its flexibility, allowing us to test different physical 
models. Moreover, an exploration of the parameter space can be 
performed in different manners and can lead to an estimate of 
the sensitivity of the fit to the different model parameters, i.e. a 
realistic error estimate. 

We applied these techniques to the special astrophysical case 
of B[e] star circumstellar environments by generating artificial 
data in order to analyse beforehand what constraints can be ob- 
tained on each parameter of the particular disc model in this 
work. The techniques will then be applied to real interferometric 
data of a sgB[e] CSE in a sequel to this paper. 

We showed in our analysis that the "geometrical" parame- 
ters such as Ri n , PAd and i can be determined with an accuracy 
£ 15 %. Mid-IR interferometric data give access to a mean tem- 
perature gradient: the temperature structure {T m and y) can be 
very well determined (within ^ 20 % and % 10 % respectively). 
It is possible to have access to the central source emission (with 
an accuracy ^ 30 %) when it has a significant contribution to the 
total flux of the object (a few % are sufficient). The remaining 
parameters of our disc model, namely n m , m and A2 are not very 
well constrained by MIDI data alone. n m is at best determined 
with an accuracy of about ^ 50 % in some cases. If A2 can be 
estimated through spectroscopic observations, then the picture 
about the n m and m determination improves somewhat. 

FRACS can be used mainly for two purposes. First, it can be 
used by itself to try and determine physical quantities of the cir- 
cumstellar matter. Admittedly, it is not a self-consistent model, 
i.e. the radiative transfer is not solved because the temperature 
structure is parameterised. From the usual habits in the interpre- 
tation of interferometric data it is nevertheless a step beyond the 
commonly use of toy models or very simple analytical models. 
This approach has indeed been very successful in the millimetric 
wavelength range (e.g. see Guilloteau & Dutrey 1998). Second, 
it can be viewed as a mean to prepare the work of data fitting 
with a more elaborate model (such as a Monte Carlo radiative 
transfer code for instance) and to provide a good starting point. 

FRACS is a tool that can help in the process of inter- 
preting and/or preparing observations with second-generation 
VLTI instruments such as the Multi-AperTure mid-Infrared 
Spectroscopic Experiment (MATISSE) project (Lopez et al. 
2006). In this respect, FRACS is not restricted to the mid-IR, and 
sub-millimeter interferometric data obtained with the Atacama 
Large Millimeter Array (ALMA) for instance can be tackled. 
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Fig. 4. Visibilities of the artificial sgB[e] circumstellar environment (model a). The visibility variations with the wavelength are shown for each 
baseline specified by the value of the projected baseline and the position angle on the sky. The circles represent the simulated observations, and 
the solid curves represent the best-fit model. 
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Fig. 5. Visibilities of the artificial sgB[e] circumstellar environment (model b). 




Fig. 6. Disc images at lOyum. (a) Image computed with the help of the Monte Carlo radiative transfer code, (b) Image of the best-fitting model 
with two power-laws (parameters of the fourth column in Table. 5) obtained with FRACS. 
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Fig. 8. Evolution of m and PAj with the inclination i. Left: maps for the couple m and PAj ; right: x 1 maps for the couple m and i. Contours 
are drawn for^f min + A^, with = 0.3, 1, 3. From top to bottom the inclination i takes the value 20 °, 50 ° and 90 °. The results correspond to 
model 4, 6, and 10. The limits of the maps have been set to ±25 % of the true values of the parameters. 



G. Niccolini et al.: Fast ray-tracing algorithm for circumstellar structures (FRACS) 



17 





- -2.9 

6 




0.0050.0100.0150.0200.0250.030 0.05 0.10 0.15 0.20 0.25 0.30 

n in (1/m 3 ) n in (1/m 3 ) 



Fig. 9. maps for the parameters «i„, Ai and in. The results presented here are those of model (a) (t = 0. 1, left part) and model (b) (t = 1, right 
part). Contours are drawn for^ min + A^ 2 , with A^ 2 = 0.3, 1, 3. The three possible maps corresponding to the combination of these parameters 
are represented. These three parameters are better constrained in model (a). 



